/*  slfdivkn.cpp  by K.Tsuru  */
// function ID = 215   DRADIX, BRADIX
/********************************************************************
SLong and SInteger classes
This is a protected member function which user cannot access.
It provides the division and remainder.
 q = m/n, r = m - q*n;    m and n both positive
It is assumed that the conditions such as m > n > 0, n > 1 have been
checked in the LLDiv(). The Knuth's method is used.
If you made a new program please test by taking
    a = 10^15-1 = 99999 99999 99999 (generally R^n-1 with radix R)
    b = 10^10-1 = 99999 99999       (id.)
    r=a*b;  c = r/a;
and check that b == c. The division "r/a" causes an "add back".
**********************************************************************/
#ifndef SN_H
#include "sn.h"
#endif
static const char* func = "KnuthLLDiv";
Ldiv_t SLong::KnuthLLDiv(const SLong& m, const SLong& n, bool needRem){
    if( (m.Head() < n.Head()) || (m.Sign() <= 0) || (n.Sign() <= 0) ){
        m.SetError(m.SYNTAX_ERR, func, 215);
    }
    int j, uh, vh;
    fType d, qq, v1, v2, rdx = m.Radix();
    SLong q(m.Type(), m.Head()+2-n.Head()), u(m);   //work area for quotient and remainder
    SLong v(n),  b1(m.Type(), 0);

//normalization and initialization
    d = rdx/( v[v.Head()] +1 );

    u = LsMult(u, d);   v = LsMult(v,d);    //u*=d, v*=d

    j = uh = u.Head() + 1;
    u.Reserve(j+1);         //for calling u(j), u(j) = 0

    vh = v.Head();
    v1 = v.figure(vh);      v2 = v.figure(vh-1);    // v1 >= rdx/2

    v.ShiftArray(uh - vh-1);    // v = v*rdx^(uh-vh-1)

    long w, w1;
    while( j > vh ){
    //It gets a tentative quotient "qq" in one figure by dividing v into u.
        w = (long)u.figure(j)*(long)rdx + (long)u.figure(j-1);
        qq = (u.figure(j) == v1 ) ? (rdx -1) : fType(w/v1); //error of qq <= 2

        w1 = ( w - (long)qq*(long)v1 )*(long)rdx + (long)u.figure(j-2);
        while( (long)v2*(long)qq > w1 ){
            qq--;
            w1 = ( w - (long)qq*(long)v1 )*(long)rdx + (long)u.figure(j-2);
        }
        u = u - LsMult(v, qq);      // u = u - qq*v;

        int us = u.Sign(215);
        if(us < 0){ //[Add back] u < 0, The probability is less than 2/rdx
// fprintf(stderr, "\"add back\" occurs.\n");
            b1.SetLong(1);
            b1.ShiftArray(j+1);     // b1 = rdx^(j+1)
            u += b1;                // a complement of b
            qq--;
            u += v;
            u.figure[j+1]=0;        //neglect the carry
            u.CheckArray(215);
        }
        q.figure[j - vh -1] = qq;
        if(us == 0) break; // added since ver 2.181. Thanks Y.Senda for the suggestion on Dec 5, 2009.
        v.ShiftArray(-1);           // v /= rdx;
        j--;
    }

    q.SetSign(1);
//It gets the figure positions.
    j = uh - vh-1;
    while(!q.figure(j)) j--;
    q.aHead = j;
    j = 0;
    while(!q.figure(j)) j++;
    q.aTail = j;

    q.CutDown(q.POP);   //need POP
    q.DoCutDown();      //reduce the size if possible

    if(!needRem) return q;  // Ldiv_t(q, 0.0);
    if( u.Sign(215) ){  //remainder!=0
        u = LsDiv(u, (ulong)d);     u.SetSign(1);   //u /= d
    }

    if( u.Verify() ) {// added since ver 2.181.
        SLong v = m - n*q -u;
        if( v.Sign(215) != 0) v.SetError(v.FATAL, func, 215);
    }
    return Ldiv_t(q, u);
}